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The Nikolaevskiy equation has been proposed as a model for seismic waves, electroconvection 
and weak turbulence; we show that it can also be used to model transverse instabilities of fronts. 
This equation possesses a large-scale "Goldstone" mode that significantly influences the stability of 
spatially periodic steady solutions; indeed, all such solutions are unstable at onset, and the equation 
exhibits so-called soft-mode turbulence. In many applications, a weak damping of this neutral mode 
will be present, and we study the influence of this damping on solutions to the Nikolaevskiy equation. 
We examine the transition to the usual Eckhaus instability as the damping of the large-scale mode is 
increased, through numerical calculation and weakly nonlinear analysis. The latter is accomplished 
using asymptotically consistent systems of coupled amplitude equations. We find that there is a 
critical value of the damping below which (for a given value of the supercriticality parameter) all 
periodic steady states are unstable. The last solutions to lose stability lie in a cusp close to the 
left-hand side of the marginal stability curve. 



I. INTRODUCTION 

The Nikolaevskiy equation 

d ^ +u ^ = I ru ~ ( 1 + 75^2 ) u I (!) 

has been widely studied, because of its application to several physical systems and its interesting nonlinear dynamics. 
It can also be written in the alternative form 





dt 2 V dx 




(2) 



after introducing u — <p x . 

The equation ([TJ was first derived, in an extended form, by Nikolaevskiy as a model for seismic waves in the earth's 
crust (TJ [2] - More recently, other applications of (JTJ or ^ have been proposed. Fujisaka and Yamada [3] and Tanaka 
[4] have derived ^ as a possible phase equation arising in reaction-diffusion systems. Turbulence in electroconvection 
was studied experimentally by Hidaka et al. [5], who proposed that the chaotic dynamics was of a type described 
qualitatively by (JT|). More generally, the Nikolaevskiy equation can be regarded as a canonical model for a pattern- 
forming system with the additional feature of symmetry under the transformation cf> — > <p + constant, in the form ^ 
[5J, or with Galilean symmetry, in the form ([!]) [TJ. Therefore, (JT[ ) is a suitable model for pattern- forming systems with 
these additional symmetries. Finally, as we show in Section [HJbclow, Q can be used to describe finite-wavelength 
instabilities of travelling fronts. 

The Nikolaevskiy equation exemplifies what is known as "soft-mode turbulence" . This form of chaotic dynamics 
arises from the interaction between a pattern of finite wavenumber, which appears for r > 0, and a long- wave neutral 
(or "Goldstone") mode. The neutral mode arises as a direct consequence of the additional symmetries discussed 
above. In a sufficiently large domain, all spatially periodic steady states are unstable, as was shown by Tribelsky and 
Velarde [8], and numerical investigations of (JT|) show irregular chaotic patterns 0171 [9]. In Ref. [7J we showed that 
in this turbulent regime, (JlJ exhibits an unusual scaling, with the amplitude of u proportional to the 3/4 power of 
the supercriticality parameter r, whereas it is more usually proportional to the 1/2 power in other pattern- forming 
systems. Fujisaka et al. [TD] extended ([2| to two spatial dimensions and found that this same scaling law holds. 

In this paper, we consider the effect of adding a weak damping term to the neutral mode in Q. The motivation for 
this is that in real experimental systems, the symmetry that gives rise to the neutral mode will often be weakly broken. 
For example, true Galilean symmetry is broken by the existence of distant boundaries. In the case of electroconvection, 
a weak magnetic field damps the neutral mode |llj . A second motivation for including the damping term is to improve 
our understanding of the appearance of chaotic dynamics in the undamped equation (JT|) . If the neutral mode is strongly 
damped, then the dynamics of patterns is governed by the Ginzburg-Landau equation and stable, steady patterns 
are observed. As the damping is reduced, a transition to Nikolaevskiy chaos can be observed and investigated. Our 
work extends the recent study by Tribelsky [12] , 
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In Sections [TT] and III we motivate and discuss the undamped and damped versions of the Nikolaevskiy equation. 
Then in Section IV we compute numerically the steady spatially periodic "roll" solutions of the damped equation 
and compute their secondary stability. This numerical calculation shows that there is an additional pocket of stable 
rolls, beyond those considered by Tribelsky [IS]; indeed this pocket is significant in containing the last rolls to remain 
stable as the damping coefficient is reduced. Most of the stability boundary of the rolls corresponds to large-scale 
modulational instabilities; we then analyse these instabilities in three regimes, depending on the relative sizes of the 
supercriticality and damping parameters. In the strong-damping regime, the relevant weakly nonlinear description 
is a modified Ginzburg-Landau equation. In the two other regimes, where the damping is either moderate or weak, 
our analyses are based on three coupled amplitude equations, for the amplitude and phase of the rolls and for an 
associated large-scale field; in each case, these amplitude equations may be reduced to a single, third-order evolution 
equation for the phase. In all three regimes, we examine the secondary stability of rolls, and in certain cases we are 
also able to determine the effects of the leading nonlinear term and comment on the direction of the bifurcation. We 
demonstrate the subcritical onset of instability in some cases, including the no-damping case, consistent with the 
observed sudden onset of the instability in numerical simulations of . 



II. THE NIKOLAEVSKIY EQUATION AS A MODEL FOR TRANSVERSE INSTABILITY OF FRONTS 

In this section we explain how the Nikolaevskiy equation may arise as a model for finite-wavenumber transverse 
instabilities of planar fronts. Consider a planar travelling front in a homogeneous, isotropic medium, arising from, for 
example, a combustion problem or a system of reaction-diffusion equations. Suppose that the front is travelling at 
speed c in the z direction and that the x coordinate is directed along the front. Consider now small perturbations to 
the uniform front. Let the position of the front under perturbation be denoted by z — ct + (f>(x,t). Now the evolution 
equation for <p must have the property that it does depend on the value of <f>, because of translational invariance; 
the equation can depend only on x-derivatives of cj>. Furthermore, reflection symmetry in x means that the only 
permissible linear terms involve even derivatives of 0. Therefore, the linearized equation for must have the form 

d<t> d 2 cj) d 4 (f> d 6 cj) 

8t= a2 dxi +a4 dx-± +a6 dxt + ---> (3) 
where the a,j are constants. The corresponding dispersion relation is 

A = -a 2 k 2 + a 4 k 4 - a 6 k 6 H , (4) 

for the growth of modes <p ~ exp(At+ifca;). From Q it is apparent that there are two possible types of instability of the 
front. If a2 < and 04 < 0, then the planar front is unstable to a band of wavenumbers near k — (Figure [I] dashed 
line). This case has been widely studied, and leads, in the weakly nonlinear regime, to the Kuramoto-Sivashinsky 
equation [T3JQ3], which, with suitable rescalings of x and t, may be written in the form 

1 fd<t>\ 2 d 2 / ,a 2 <A 

: =-1^2 \ r ^+j^2 • ( 5 ) 



dt 2 \ dx J dx 2 \ dx 2 _ 

But if 02, 04 and ae are all positive, then long- wave modes are stable and a finite-wavenumber instability is possible, 
as illustrated by the solid line in Figure [I] After a rescaling of x and t, ^ then corresponds to the linear terms in 
([2]) . These two possible types of instability are sometimes referred to as "type II" and "type I" , respectively [23 [TB] . 

In either case, a nonlinear term must be added to the equation to limit the growth of the instability; its derivation 
is the same for either type of instability, and, as with the linear terms, only x-derivatives of </> may appear. The 
nonlinear term can be derived by a simple geometrical argument, considering the propagation of a tilted front |13j . 
Such an argument leads to a nonlinear term which in our scaling takes the form (d(j)/dx) 2 /2, as in ^ and 

We have shown that ([2| may be regarded as a model equation for the transverse instability of a planar front at finite 
wavenumber. More generally, it is apparent that any physical system that gives rise to the Kuramoto-Sivashinsky 
equation might also lead to the Nikolaevskiy equation, depending on which of the two possible forms shown in Figure[T] 
is taken by the dispersion relation. Another example [31 HI HZ] is the derivation of nonlinear phase equations; this leads 
to the Nikolaevskiy and Kuramoto-Sivashinsky equations, in different parameter regimes, as the dispersion relation 
changes from type I to type II. The argument in the case of phase equations is essentially identical to that given above 
for fronts, following from the invariance under addition of a constant to the phase, and under reflection in x. 

There are many studies of the transverse linear stability of fronts in the literature, arising from a variety of different 
physical and chemical systems. A review of the literature reveals several examples of a finite-wavenumber (type 
I) instability, consistent with the Nikolaevskiy equation. For example, in their study of an evaporation front for 
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FIG. 1: Curves illustrating two possible dispersion relations for front instabilities. Solid line: finite- wavenumber instability, 
leading to Nikolaevskiy equation (type I). Dashed line: long- wave instability, leading to Kuramoto-Sivashinsky equation (type 
II). 

condensed matter heated by a laser, Anisimov et al. |18j found a linear spectrum of type I for transverse instability. 
Both types of instability were found in an experimental and theoretical study of a reaction front between two chemicals 
in a Hele-Shaw cell [EJ QH [20] . 

III. THE DAMPED NIKOLAEVSKIY EQUATION 

The complex chaotic behavior of the Nikolaevskiy equation at onset arises from the interaction between the unstable 
modes near k = 1 and the neutral mode at k = 0. However, in real physical systems it is likely that the symmetry 
giving rise to the neutral mode (such as translation or Galilean symmetry) will be weakly broken. For example, in 
the case of a propagating front, considered in Section [TTJ the translation symmetry may be broken by the presence of 
boundaries, or by a slight variation in the basic state with position. With such a weak symmetry breaking, the k = 
mode will not be truly neutral, but will be weakly damped instead. 

A further motivation for including a damping term in the Nikolaevskiy equation is that it allows us to study 
the transition from ordered stable patterns to soft turbulence. Such a probing of the origins of the turbulent state 
is not possible in the original Nikolaevskiy equation (JlJ, since chaos is observed directly at onset (provided the 
domain is sufficiently large). Introduction of a damping term gives another parameter, which can be used to control 
and investigate the onset of chaos. A similar approach has been used to examine the dynamics of the Kuramoto- 
Sivashinsky equation [5TJ [55] . 

We thus consider in this paper the "damped Nikolaevskiy" equation 

du Ou d 2 ( ( d 2 \ 2 \ / 3 2 \ 2 

dt +U d X = -d X 2 { rU -{ 1+ d X 2 ) M J-T + ^J U > (6) 

with damping coefficient v > 0; in particular we shall be concerned with the dynamics of ^ near the onset of pattern 
formation. In examining the stability of roll solutions, we shall find it useful to consider various cases for the relative 
sizes of the supercriticality parameter 

r = e 2 

and the damping coefficient v, and hence we write 

v = e>, (7) 

where fi = 0(1); we shall examine below a variety of pertinent values for s. 

Of course there are many ways in which one might add damping to the Nikolaevskiy equation. The damping term 
in ([6| is chosen so that all modes are linearly damped (with the exception of the mode with wavenumber k = 1); 
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large-scale modes decay at a rate —v. Furthermore, significant analytical simplifications follow from the fact that in 
([6]) the onset of linear instability is independent of v (i.e., the critical value of r is r c = and the critical wavenumber 
is k c = 1, for any v > 0). Note that our choice of damping term differs from that of Tribelsky [IS], who considers 
instead a damping term proportional to — u — d 2 u/dx 2 . Consequently, in his formulation both r c and k c become 
functions of v, leading to additional algebraic complications. However, the two formulations are equivalent under 
appropriate rescalings of x and t, and appropriate mappings between the two parameter sets. 



IV. SECONDARY STABILITY OF ROLLS: NUMERICAL RESULTS 

We begin our investigation of ^ by numerically computing spatially periodic steady "roll" solutions, and deter- 
mining their secondary stability. To do so, first we fix e and compute the roll solution u(x) for a given choice of k and 
v by integrating (|6) forwards in time in a box of length 2n/k. In our pseudospectral numerical code to achieve this, 
u is approximated by the truncated Fourier series 



JV 

„mkx 



-N 

Disturbances u'(x) to this roll solution then satisfy, from ([6]), 

du' .du du' d 2 ( , / d 2 \ 2 A / d 2 \ 2 , 

-m +u d x + u ^ = -d^[ ru -[ 1 + a^) u )-»[ 1 + dx- 2 ) < (8) 

and these disturbances are sought in the form 

N 

u'(x) =e at+lpx ^u' n e nlkx (9) 
— JV 

(where the sum is again truncated for numerical purposes). Equation Q is thus reduced to a matrix eigenvalue 
problem for the growth rate a. By computing a{p;k 1 t,v) for p in the range —k/2 < p < k/2, we are then able to 
determine the stability of the roll solution u(x). 

Results for e = 0.1 and e = 0.05 are summarized in Figure [2] where rolls are stable inside the regions bounded by 
solid curves. For large enough damping (i.e., towards the top of each plot), the familiar Eckhaus stability boundaries 
are recovered. For sufficiently small damping, by contrast, all rolls are unstable - this shows that the known instability 
of all rolls in the undamped Nikolaevskiy equation extends to the case of (sufficiently small) finite damping, for a fixed 
value of the supercriticality parameter. Intermediate numerical results show the transition between the two limits. 

Several features of Figure [2] merit comment. First, the region of stable rolls, which is roughly symmetrical about 
k = 1 at large damping, becomes distinctly asymmetrical as the damping is reduced; rolls with k > 1 tend to be less 
stable than those with k < 1. Second, stability is finally lost, as v is decreased, in a small cusp near the left-hand side 
of the marginal stability curve. For e = 0.1, this corner is at v w 0.02 and k « 0.950 (the roll existence boundary is, for 
this value of u, at k « 0.9493); for e = 0.05, we have found it much harder to compute the very thin cusp right down 
to its tip. The left and right sides of the cusp correspond, respectively, to monotonic and oscillatory long-wavelength 
instabilities. In fact the entire boundary corresponds to long-wavelength instabilities (|p| <C 1), with the exception of 
the segment marked "/", where the instability has finite wavelength (p — O(l)). For even smaller values of e than 



shown in the figure, this finite- wave length segment is absent, a point to which we shall return in Section VC 1 

A more conventional picture of the stability balloon for rolls is shown in Figure [3] where the damping parameter 
v is fixed to be 0.02 and the stability of rolls is shown in the (fc,r)-plane. Again rolls are stable inside the solid 
curves. The region of stable rolls lies almost entirely in the region k < 1. For r > 0.003, all rolls are unstable, except 
those exceedingly close to the left-hand part of the marginal stability curve; these rolls are last to lose stability as r 
is increased. Finally, for r > 0.01, all rolls are unstable. A corresponding picture for smaller damping coefficient v 
would show stable rolls in a similar, but smaller region, confined closer to the point of onset (r = 0, k = 1). In the 
next section, we shall argue on the basis of asymptotic arguments that in the limit v — > (approaching the undamped 
Nikolaevskiy equation), the size of the stable region shrinks to zero. Crucially, however, we expect a small, but finite, 
region of stable rolls in the (fc, r)-plane for any nonzero value of the damping coefficient v. It is only for the undamped 
case v = that all rolls are unstable at onset. 
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FIG. 2: Stability boundaries of rolls in q6b 9 computed numerically for (a) e = 0.1 and (b) e = 0.05. Rolls with wavenumber k 
exist between the dashed curves, and are stable inside the region enclosed by solid curves. The annotations on the stability 
boundary indicate the type of instability suffered by the rolls upon crossing that part of the boundary: modulational steady 
("s"), modulational Hopf ("ff") or finite- wavelength ("/") instability. The crosses correspond to numerical simulations of the 
initial-value problem |6| in a large computational box and indicate the most extreme values of the wavenumber for which rolls 
are found to be stable to small perturbations (necessarily restricted to those that fit into the computational box); that all 
crosses lie inside the roll stability boundary provides a consistency check on our results. 
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FIG. 3: Stability boundaries of rolls in computed numerically for v = 0.02. Rolls with wavenumber k exist between the 
dashed curves, and are stable inside the region enclosed by solid curves. 

V. SECONDARY STABILITY OF ROLLS: ANALYTICAL RESULTS 

The numerical results in the previous section reveal a number of different instabilities of rolls in the damped 
Nikolaevskiy equation. In this section, we consider in detail three scalings for v in (JsJ) which shed particular light on 
the numerical stability results illustrated above: these correspond to choosing in turn s = 1, s = 3/2 and s = 2 in d7j. 

A. Strong damping: s = 1 



The case s = 1 represents strong damping, in the sense that the decay rate v = e/i of the large-scale mode is O(e), 
which is much greater than the 0(e 2 ) growth rate of the pattern-forming mode. The effect of this strong damping is 
simply to modify the usual Ginzburg-Landau amplitude equation and thus modify the Eckhaus stability boundary, 
as we now demonstrate. Our conclusions below are consistent with those obtained by Tribelsky |12j . although his 
analysis proceeds directly from a substitution of |9]) into ^ , rather than the amplitude-equation framework developed 
below. 

We apply to (T6j) the usual weakly nonlinear scalings and expand the solution as 



e{A(X, T)e ix + c.c.) + e 2 u 2 + e 3 u 3 + 



(10) 



where X — ex and T = e 2 t, and where c.c. denotes the complex conjugate of the preceding term. At O(e), ^ is 
satisfied. At 0(e 2 ), we find 

u 2 = (-JLi A 2 : e 2lx + c.c.) + f(X,T), 

where the large-scale mode f(X,T) is arbitrary. At 0(e 3 ), an equation arises for the mean mode /, which represents 
a balance between the driving term (|^4| 2 )x and the damping term —fif: 



0=-(\A\ 2 ) x - f if. 



(11) 



Note that the linear terms fx and fxx do not appear at this order. The governing equation for A is also obtained at 
0(e 3 ), by applying the usual solvability condition. It turns out to be 



A T = A + AA] 



ijcx- ±\A\*A-ifA, 

from which / may be eliminated using (11) to give a single amplitude equation for A in the form 

A T = A + AA XX - ±\A\ 2 A + i(\A\ 2 ) x A/^. 



(12) 
(13) 
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FIG. 4: Stability boundaries of rolls according to (191, in (q, ^t)-space, where the damping coefficient v 



e/j. Rolls exist for 
0, the stable rolls 



— 1/2 < q < 1/2 and are stable in the region between the two curves, where indicated. In the limit /i 
are those with q < 0. In the limit as fi — > oo, the usual Eckhaus stability boundary is recovered, so that rolls are stable for 
-(1/12) 1 / 2 < q < (1/12) 1 / 2 . 



This is the usual Ginzburg-Landau equation for A, but with an additional nonlinear derivative term arising from the 
damping. Extensions to the Ginzburg-Landau equation including terms such as i(|A| 2 )x^4 are well known [23] [24] . 
but the different scalings involved here have the consequence that no further terms need be included to ensure that 
all terms of a given asymptotic order are present: with our scalings, (13 1 is complete. Note that in the limit fi — > oo, 
we recover the usual Ginzburg-Landau equation for A. 



It is straightforward to analyse the stability of patterns in (13 1 - see Mancebo and Vega [25]. Steady patterns with 
A = a a expiqX exist for q 2 < q 2 = 1/4, with the real amplitude do = 6(1 — Aq 2 ) 1 / 2 . The stability of these solutions 
may be determined by writing A — (ao + a(X, T)) expiqX and linearizing in the perturbation a, to give 



a T = -^a 2 l (a + a*) + 8iqa x + 4a X x + ial(a x + a* x )/n. 



After writing a = b + ic, and separating ( 14 1 into real and imaginary parts, we find that b and c obey 



^-n 2 b - 
18 « o 



,2i 



bx '- 

c T = Sqb x + 4:C X x + 2af ) b x /^- 
Then by supposing b and c each to be proportional to exp(AT + llX), we find that the growth rate A satisfies 

A 2 + (iog + 8Z 2 )A + 16Z 4 + §a§Z 2 - Uq 2 l 2 - 16ql 2 al/fi = 0; 
only monotonic instability is possible, and there can be instability if and only if 

l 2 (72l 2 + a 2 - 288q 2 - 72qa 2 ,/ » < 0. 



It is clear from (18 1 that the most "dangerous" modes are those with small I, so the pattern is stable if 

12g 2 ) > 72g(l - Aq 2 )/^. 



(1 

The corresponding stability boundaries are shown in Figure [4] 

In the limit fi — + oo, the usual Eckhaus stability condition (that rolls are stable for q 2 



< q, 



2 - io 2 
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(14) 

(15) 
(16) 

(17) 

(18) 

(19) 
) is 



recovered from (19). As \i is decreased, however, an asymmetry develops in the stability region: patterns are stable 



for f(fi) < q < g{fJ>), where f(fi) € [— h, —q e ) and g(n) € [0,q e ) are each increasing functions of /i. When /i is small, 
the stability condition is roughly g(l — Aq 2 ) < 0, so that patterns with q > are unstable and patterns with q < 
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FIG. 5: Space-time plots of numerical solutions of ( 13 1, for /j, — 0.1 (left, real part of A plotted), and the undamped Nikolaevskiy 
equation Q with r = 0.01 (right). 

are stable. A corresponding conclusion is reached by Tribelsky (whose results are plotted in his Figure 1(a)). If we 
compare the results of this analysis with the numerical stability calculations described above in Section |IV| we see 
that the stability boundaries in Figure [4] correspond, roughly, to those in Figure [2] above v w 0.2 (for e = 0.1) and 
above v w 0.07 (for e = 0.05). 

The above analysis shows that for small /i, patterns with q < are stable in (|13|. However, numerical simulations 



of (13 1 in a periodic domain, started from a random, small amplitude initial condition do not always lead to stable 
patterns. A typical simulation, for /i = 0.1, in a domain of size 100, is shown in Figure [5] (left). The solution is 
dominated by a state with 8 rolls in the domain, with a wavenumbcr q — —8 x 27r/100 » —0.5026, which is just 
outside the range of wavenumbers for which steady patterns exist. This mode must therefore decay, until a bursting 
event occurs, involving other wavenumbers, from which the 8-roll state again emerges. Simulations with different 



domain sizes show that (13) seems to exhibit a preference for modes near to the left-hand limit of the region of 
existence of steady states, q = —1/2. Note that this behavior is shared by the damped Nikolaevskiy equation pj, 
as shown in Figure [2] Another qualitative similarity between the behavior of (13 1 and the Nikolaevskiy equation is 
that both exhibit irregular bursting behavior (albeit with rather different detailed structures). For comparison, a 
simulation of the undamped Nikolaevskiy equation Q in a domain of size 200, for r = 0.01, is shown in Figure [5] 
(right). 



In the limit fj, — > 0, it is clear from (17 1 that the eigenvalues A become large, of order pT x l 2 . This indicates that 
a different scaling is required in which the damping v is smaller than O(e) and the growth rate of the instability is 
greater than 0(e 2 ); this next scaling is considered in the following section. 

B. Intermediate damping: s — 3/2 

The analysis above corresponds to damping that is sufficiently strong to enslave the large-scale mode to gradients 
of the pattern amplitude. More subtle effects of damping can be explored if it is rather weaker: if s = 3/2 in Our 
starting point is now the general roll solution of ^ with v — e 3 / 2 ^, which may be written in the form 

u = ea e {1+tq)vx + c.c. + 0{e 2 ), 

where a = 6(1 - V) 1 / 2 . 

A proper treatment of the stability of rolls requires consideration of the evolution of both amplitude and phase 
of the rolls, together with a large-scale mode. All three modes couple together, and their relative scalings, together 
with appropriate length and time scales for their evolution, are chosen below to enable a balance in the linearized 
perturbation problem; this is essentially the balance considered in [TJ. However, by making appropriate absolute 
scalings for the three perturbation quantities, we are able to extend the linear results in [7J and accommodate the 
first nonlinear term in the evolution equations that follow. We shall see below that this enables us to describe the 
sub- or supercritical nature of the onset of instability of the rolls. It turns out that the appropriate scaling for the 
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three perturbation quantities is accomplished by writing 

u = e(a + e 1 / 2 o(X, T ))e lx e lqex e ltl/4 ^ X < T) + c.c. + • • • + e 7/4 /(X, T) H , (20) 

where now 

X = e 3 / 4 x, T = e 3 / 2 i. 

Note that in order to simplify notation we use the same symbols (X and T) in different sections to represent different 
length and time scales. Notation within an y on e section is consistent, so no confusion should arise. 

Then after a systematic substitution of (20 1 in ^ and a consideration of terms at successive powers of e 1 / 4 we 
eventually find 



dT 
d£ 
dT 
()a 
df 



dX 2 
d 2 f 
OX 2 
, d 2 a 



dX 2 



-f, 
- 4a n 



2a 



da 
dX' 



dX 



dX' 



(21) 
(22) 
(23) 



The linear terms in this equation may be identified with the linear system considered in [7J , with the addition here of 
the new damping term —)J,f in ( |22[ ) (the quantities b, c and / in [7] correspond, respectively, to our a, ao(j) and /); the 
nonlinear term was not computed in [7 . These three equations may then be reduced to the nonlinear phase equation 



d_ 
df 



d 2 



dX 2 



d_ 
df 



d 2 



dX 2 



-IQai 



dcj)\ d 2 (j> 
q+ dX dX 2 ' 



(24) 



In analysing (24) 
-16alq(/)xx 



it is useful to consider first the linearized problem, for which the right-hand side is si mpl y 
Th is term re present s the coupling between <j>, a and /: in its absence, no instability can result in (24 1. 
Thus in Sections VB 1 and | V B 2 below we assume a^q — 0(1) 
below must be reconsidered when \a 2 1 q\ is small (cf. [8]) 



However, we should note that the analysis described 
we see immediately that there are two circumstances in 
which a separate consideration is warranted (a <C 1 and \q\ < 1); we shall return to this important point below, in 
Section EC] 



No damping 



We begin by summarizing relevant results in the absence of damping. When /i = we recover from ( 24 1 the 
secondary stability results obtained by Tribelsky and Velarde j8] and elsewhere by ourselves [Jj: there is monotonic 
instability of the rolls with q > 0, and such rolls are unstable to disturbances with wavenumbers I in the range 
< \l\ < l cm = (a^) 1 / 4 ; there is oscillatory instability for q < 0, and here the instability strikes for < \l\ < l co = 
(-2a^/25) 1 / 4 . Thus all rolls are unstable at onset in the undamped problem. 

It is instructive to extend the linearized analysis by considering the weakly nonlinear development of these insta- 



bilities, according to (24 1. We shall focus on computing the coefficients of the nonlinear terms in Landau equations 
for the amplitudes of the disturbances; the signs of the real parts of these coefficients will indicate the direction of 
bifurcation to the perturbed state. We consider first the monotonic instability for q > 0, and examine the evolution 
of a disturbance whose wavenumber satisfies \l — l crn \ = 0{8), where < <5 <C 1. Expanding 

cb = 6fa + 8 2 <j> 2 + • • • , (25) 



where 



and T 2 = S 2 T, we find at 0(5 3 ) that 



A(T 2 )sinlX 



dT 2 180g ' 
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FIG. 6: Stability boundaries of rolls, based on analysis of the case of intermediate damping, v = e 3//2 /i. Rolls exist between 
the dashed lines and are stable or unstable as indicated. The solid vertical straight line represents the onset of a monotonic 
bifurcation and the solid curve represent a Hopf bifurcation. 
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Thus the bifurcation is always subcritical and hence leads to the observed explosive growth of the instability in 
simulations of the undamped Nikolaevskiy equation [71 [5]. For the oscillatory instability for q < 0, we instead suppose 
that \l — l co \ = 0(S); again expanding <fi as in (p5j), but now with 



01 = c lujT (B(T 2 )e llx + C(T 2 )c- llx ) + c.c. 
where w — (— 48a 2 ,<7/25) 1 / 2 , we find at 0(5 3 ) the Landau equations 



—q dB 
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alBdT 2 


I 72639 


8071 


-q dC 
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a 2 CdT 2 
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(27) 

Since q < 0, both travelling waves and standing waves branch supercritically, with the travelling wave stable. In 
practice we expect the explosive development of the monotonic instability to dominate the smooth onset of the 
oscillatory instability in simulations of the undamped Nikolaevskiy equation, except under extremely controlled (and 
contrived) conditions. 



2. Nonzero damping 



For /i > 0, we again begin by examining the linearized behavior of (24 1. A monotonic instability again occurs for 
all rolls with q > 0, with the stability margin being given by l\ m + [il\ m — qa^ = 0. It thus follows that 

'<UO) = Q a - V l lm < Q a = 'cm(°)> 

and hence the rolls are unstable to a smaller range of disturbance wavenumbers than in the undamped case. An 
oscillatory instability arises for rolls with q < qo, where qo = — /i 2 /(2a§); such rolls are unstable to disturbances with 
< |Z| < l COl where now 5Z 2 G = —fj, + (— 2ga 2 ,) 1 / 2 . We may rewrite the condition q < qo for instability as 

72g(l -Aq 2 ) < -fi 2 . (28) 



Then since the cubic function in (28) has a minimum value of —8y/3 at q = — V3/6, this oscillatory instability can 



occur only for /i 2 < 8v3, that is, for sufficiently weak damping (this is consistent with no oscillatory instability being 
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observed in the previous scaling of stronger damping in Section V A ) . In the limit of small /z the stability boundaries 



in p8| approach q — and q = —1/2, so that all rolls with q < are predicted to be unstable to an oscillatory 
disturbance, just as in the undamped case. 

The stability boundaries for this scaling are illustrated in Figure [6] Note that the curved line of the Hopf bifurcation 
is consistent with the numerical stability results for the full damped Nikolaevskiy equation (|6| shown in Figure [2] 
(except, of course, that the analysis above does not capture the finite- wavelength instability denoted by "/" in 
Figure [2]). Figure [6] suggests that two regions of stable rolls, near q — and q — —1/2, extend right down to zero 
damping. But in each of these cusped regions, \a 2 ] q\ = |6q(l — 4g 2 ) 1 / 2 | <C 1, and we recall that the present scaling is 
not valid in this limit. Thus we should disregard the two cusps in Figure [6] at small damping and instead examine more 
closely the stability of rolls in this limit; to resolve the correct behavior of the stability boundary at small damping, 
we consider in the next section a further (and final) scaling for v. 



C. Weak damping: s = 2 

The final small-e asymptotic description of the roll stability boundary may be described by setting s — 2, so that 
v = e fi. In this case, the damping rate of the mean mode is of the same order as the growth rate of the pattern- 
forming mode. There are two cases to consider where the scaling of the previous section breaks down, corresponding 
to \q\ <C 1 and a <C 1. 



1. Central corner (\q\ <C 1) 



As Tribelsky and Velarde [8] have pointed out for the case of no damping, the scalings of Section VB break down 
when q is small. To resolve the small-g behavior, we adopt Tribelsky and Velarde's scalings (but in an amplitude- 
equation framework) and consider 

u ~ e(a + e 2 a(X, T))e lx e le ' ' qx e I£< ^ < T) + c.c. + • • ■ + e 3 f(X, T) + ■ ■ • , (29) 



where now X — ex and T = e t and ao — 6. Note that whereas (20 1 concerns basic roll wavenumbers k = 1 + 0(e), 
(129 1, by contrast, concerns the much narrow wavenumber band k = 1 + 0(e 2 ). We find, after a substitution of (29 1 
in |6|, 

«*=4**-/ (30) 
dT dX 2 J ' ( ' 

df d 2 f da 

df = dx 2 -^- 2ao dx' (31) 

da d 2 a 1 2 /ll . „ 9 2 \ ^ / <90\ 2 df 

a a + a [ —a - 8q + 12— — \ — -4a \ — ) -a — . (32) 



dT dX 2 18 u U V54° y dX 2 J dX u \dX J u dX 
These three equations may then be reduced to the single nonlinear phase equation 

d_ d 2 \ ( d d 2 1 2 \ / d (f_ \ 

df ~ dx 2 J \df dx 2 + i8 a °) \df~dx 2+fi ) (l) ~ 

2a H 54 - 89 + dr + "dx 2 )a£ ~ 16a °af 3A§- (33) 



We now consider the linearized version of ( 33 ) and examine the dispersion relation for infinitesimal disturbances 
proportional to exp(aT + 1IX). It is helpful first to recall the stability results of Tribelsky and Velarde j8] for the case 
/1 = 0, which are summarized in Figure [7j and which are a special case of the more general analysis to be presented 
below. As can be seen in the figure, all rolls are unstable for /1 = 0. 

For general values of fi, rolls undergo a monotonic instability (<r = 0) when q = q m (Z,/x), where 

72q m (l 1 n) = 2l 4 + (2/x - 7l)l 2 + 66 + fi. 

Since <7 m (i,/u) > q m (l,0), it follows that the stability margin tends to move to the right as the damping is increased, 
tending to stabilize the rolls (see Figure [8|. Rolls undergo an oscillatory instability (er = iuj) when 

fi(2 + n) 2 ( V 3/i 91 \ /5/x 677 \ 2 25 



288 V 72 16 144 J \ 36 288 / 72 ' y ' 



12 




FIG. 8: Stability boundaries of rolls with damping coefficient v = e 2 /i, according to (331 
right are the monotonic stability boundaries 



for values of [i as indicated. To the 
to the left are the oscillatory stability boundaries. As the damping coefficient 



increases, the rolls become less unstable. For sufficiently large fi, a stable band of wavenumbers emerges. 



with onset frequency given by 

uj 2 = 2Al A + 2(4/i + 41)Z 2 + 2/i. 

The oscillatory stability boundary is plotted in Figure [8] for various values of /i. Damping is seen to shift the 
oscillatory stability curve to the left, again stabilizing the rolls. A further effect of nonzero damping is to stabilize 
rolls to oscillatory disturbances in the limit that the perturbation wavenumber I — ► 0. 

As can be seen from Figure [8] with increasing damping the two stability boundaries separate and eventually part, 
allowing a small band of stable rolls. We find that stable rolls exist for > /i c « 8.445 (the corresponding critical 
value of v is v c — ^ c e 2 ); thus, for a given value of the supercriticality parameter, there is a threshold value of the 
damping coefficient to allow stable rolls near k = k c , as previously shown by Tribelsky [12 . Figure [9] shows the region 
of stable rolls in (/z, <7)-space (thus the apparent small-g cusp in Figure [6] which appears to extend down to zero 
damping, is in fact a wedge terminating at finite damping) . 

The results above seem at odds with the numerical results presented in Figure [2] since there is no central wedge 
in the numerical stability boundary. However, it turns out that for the values of e used in Figure [2] the small-*? 
wedge is masked by an additional finite-wavenumber instability. For a smaller value of e, the theory in this section 



does indeed correctly predict the shape of the stability boundary, as shown in Figure 10 for e = 0.01. (It becomes 
increasingly challenging to compute the entire stability boundary numerically, so only the relevant part is presented for 
this value of e.) Rolls are stable in the wedge down to v ss 8.445e 2 , and the last rolls to destabilize have wavenumber 
approximately 1 — 4.049e 2 . These results agree with corresponding results of Tribelsky [T2] , 
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FIG. 9: Stability boundary for rolls with wavenumber 1 + e q, according to (33 1, for weak damping, with v = e fx. Rolls 



are stable inside the wedge. At the right-hand boundary they become susceptible to a monotonic instability, at the left-hand 
boundary an oscillatory instability. The wedge terminates at the point (fj,, q) ~ (8.445, —4.049). 



V 



0.003 r- 
0.0025 - 

0.002 - 
0.0015 - 

0.001 - 




0.998 0.999 



FIG. 10: Stability boundary for rolls in computed numerically for e = 0.01, near the central wedge region corresponding 
to Figure [9] Crosses indicate numerical results and the solid line corresponds to our analytical results (cf. Figure |9j. Note that 
for this value of e (and presumably also for even smaller values) the numerical stability boundary contains this central wedge 
rather than a section corresponding to that marked "/" in Figure [2] 



2. Left-hand corner (ao 1) 

We now examine the second case in which the asymptotic results of Section |V B| break down: a <C 1. This regime 
proves to include the last rolls to become destabilized as the damping is reduced. We thus investigate the stability 
of rolls close to the marginal stability boundary. Recall from Section VB and in particular Figure [6] that we expect 
stable rolls to persist to small values of the damping near the left-hand side of the marginal stability curve. 

We introduce the notation k = fc~ , fc+ for the left- and right-hand marginal stability boundaries, respectively (that 
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is, for given values of r and v, rolls exist for wavenumbers k in the range fc~ < k < fc+). We shall begin by considering 
both left- and right-hand parts of the marginal stability curve, since both lead to ao <C 1. Our analysis will then show 
that only the left-hand part of the marginal curve is relevant, being capable of supporting stable rolls. We find, from 
substitution in the equation fc 2 „[e 2 — (1 — kf n ) 2 ] — e 2 /i(l — fc 2 j) 2 = for the marginal curve, that 

fc™ = l±^-ie 2 ±(i-i M )e 3 + ---- 

It turns out that the correct scaling to resolve the stability boundary of the rolls is obtained by setting \ k— k^\ = 0(e 2 ), 
in which case |u| = 0(e 3 / 2 ), with the relevant space and time scales for the evolution of perturbations being X = ex 
and T = e 2 t. We consider here only the linearized secondary stability problem and we are thus at liberty to scale 
arbitrarily the absolute sizes of the three perturbation quantities, keeping their relative sizes fixed. Thus we suppose 
that 

u ~ e 3 / 2 (a + a{X, T ))e lkx e^ x ^ + c.c. + ■ ■ • + e 2 f(X, T) + • • • , (35) 
where the wavenumber of the rolls under consideration is just inside the marginal curve, so that 

k — k m -p e k. 

A weakly nonlinear calculation of the roll solution itself shows that the amplitude and wavenumber are related by 



= 144k. 



(36) 



Then a consideration of the terms in ^ which are linear in the perturbation quantities <f>, a and / shows that their 
evolution is governed by 



no 



dT 

df_ 
dT 
On 
df 



4a r 



2 



dX 2 



a /±4 



d 2 .f 



dX : 
da 



. d 2 a d<j) 

A dx 2TAa °dx- 



These three equations may then be combined to yield a single linearized phase equation in the form 

2 



d_ 
df 



d 2 

dX 2 



d_ _ d 2 

df dX 2 



1G 



d 2 

dX 2 



2 ^ 

dX 2 ' 



For modes proportional to exp(crT + ilX), this gives the dispersion relation 

(cr + I 2 + /i) \(a + 4l 2 ) 2 - 161 2 ] = ±8a 2 l 2 . 



(37) 
(38) 
(39) 

(40) 

(41) 



In the limit ao — > 0, (41 ) gives the growth rates a = — (fi + I ) < 0, a — 4Z(1 — /) and a — 4Z(— 1 — I), and hence 



all rolls are unstable (to a monotonic disturbance) sufficiently close to the marginal curve. We gain some insight 
into the possible existence of stable rolls near the marginal curve by next considering the limit of small perturbation 
wavenumbers, I — > 0. For fixed ao and fi, in the limit I — > 0, the leading-order balance in (41 1 is 



8l 2 (2±a 2 /n). 



(42) 



Thus near fc = fc+ rolls are unstable, to monotonic disturbances, since (2 + a 2 ,/ '//) > 0. No further analysis of this case 
is necessary: there is no pocket of stable rolls possible near the right-hand marginal stability curve; this conclusion is 
consistent with Figure [2] 

Near k — fc~, by contrast, the conclusion of monotonic instability follows from (42 1 only if a 2 , < 2/i (i.e., sufficiently 



close to the marginal curve). When a^ > 2/i, (42 1 instead indicates oscillations, and we must go further in our 



consideration of the small-/ problem in order to determine the stabi lity of the rolls. By carrying out a small-/ expansion 



of the growth rate in the form a = ail + a^l 2 H , we find from (41 ) that a\ = 8(2 — a\f fi) and ct 2 = 4(aQ//i 2 — 1). 

Thus when 



2fi < al < /i 2 , 



(43) 
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rolls are stable in the small-? limit (the lower threshold for a 2 , corresponding to a monotonic instability, the upper 



threshold to an oscillatory instability). The relation (43 1 indicates that rolls enjoy this region of stability only for 



sufficiently large damping: fi > pb c = 2. Note that since 2 < 8.445 this cusp contains the last stable rolls as the 
damping is decreased. 

In terms of the original variables in the damped Nikolaevskiy equation (JsJ) , we thus predict (for small v) the cusp 
to be at v ~ 2r = 2e 2 . We may readily check these analytical results against the numerical secondary stability 
calculations of Section [IV] For e = 0.1, as in Figure [2j we thus predict the apex of the cusp to lie at v = 0.02 (which 
we do indeed find numerically). Similarly, for e = 0.05 we predict the apex to lie at v = 0.005, and the smallest 
value of v for which we have been able to compute the cusp is v — 0.006, which provides reasonable agreement with 
the theory. (As e — > it becomes exceedingly difficult to compute accurately this part of the numerical secondary 
stability boundary.) 



For completeness, we note that we have also considered numerically the full dispersion relation (41 1 for arbitrary 
values of I: we find that (at least for values of clq and fj, up to 4, for which we have carried out calculations) the 
stability boundary is indeed determined by the small-/ expansion. 

VI. NUMERICAL SIMULATIONS 

In this section we present full numerical simulations of the damped Nikolaevskiy equation (|6]), illustrating the 
transition from regular patterns to soft-mode turbulence as the damping is reduced. The simulations use periodic 
boundary conditions and employ a Fourier spectral method for the spatial discretization. The time-stepping is an 
explicit second-order Runge-Kutta form of the exponential time differencing method |26j that computes the stiff linear 
part of the equation exactly, permitting O(l) time steps. The size of the domain is 2007T, allowing a resolution of 0.01 
in wave number. 



In the simulations shown in Figure 11 the initial condition comprises rolls with wavenumber k — 0.96, plus a small 
random perturbation. In each case, the simulation was run for a long time to remove transient features, and then 
plotted over a subsequent relatively short time interval. Note that these time intervals are not the same in each case. 
The driving parameter was fixed at r — 0.01 and the sequence of simulations shows the effect of reducing v. For these 
parameters, all rolls that fit into the computational domain are predicted to be unstable for v < v c where v c w 0.057. 
The bifurcation at v — v c is the Hopf bifurcation whose stability boundary is shown in Figures [2] and [6] 

The simulations confirm the predicted bifurcation type, since for v < v c a small oscillatory modulation to the rolls 
grows. This bifurcation seems to be supercritical since the oscillations equilibrate to a small periodic modulation of 
the rolls. For v = 0.04 this oscillation takes the form of a standing wave, as shown in Figure [TX| (top). At v — 0.03, the 
solution is no longer periodic in time and shows occasional irregular bursts of instability. In this case, the dominant 
wavenumber is k = 0.95, at the left-hand boundary of the roll-existence region (i.e., near k — k~ n ), consistent with 



the analysis of the preceding section and the behavior of ( 13 1. Also, counter-intuitively, the amplitude of the solution 



becomes significantly lower as the damping is decreased. This is because of the transition from the 0(e) scaling in 
the strongly damped case to the 0(e 3 / 2 ) scaling [7] in the undamped, unstable regime. For v — 0.02, the bursts are 
more frequent and a broad spectrum of wavenumbers is present in the solution. Finally, the lowest panel of Figure [TT] 
shows the case of zero damping. This is significantly different from the weakly damped case, showing grain boundaries 
between regions of travelling rolls. 

VII. DISCUSSION AND CONCLUSIONS 

The Nikolaevskiy equation is an important model of a wide range of physical systems, including certain convec- 
tion problems, phase instabilities and transverse instabilities of fronts. It arises naturally in systems with a finite- 
wavenumber instability and a translation symmetry for the dependent variable. However, in many applications it is 
likely that this symmetry will be weakly broken. 

The damped Nikolaevskiy equation (|6]) allows us to investigate the effects of weak symmetry breaking, and describe 
the transition between the appearance of soft-mode turbulence at onset when v = and the more gradual development 
of complex dynamics typical of damped systems when the damping coefficient v > 0. Our numerical investigation of 
the secondary stability problem for steady roll solutions shows how the instability of all roll states in the undamped 
case evolves into the more common Eckhaus scenario, whereby rolls with wavenumbers sufficiently close to critical 
are stable, when the damping is sufficiently great. If we fix the supercriticality parameter r, then it follows from the 



asymptotic results of Section V C 2 that there is a critical value of the damping coefficient v = v c (given by v c ~ 2r 



when r is small) below which all rolls are unstable; for v > v c some rolls are stable. 




FIG. 11: Space-time plots of the numerical solution of <[6| at r = 0.01. Top to bottom: v = 0.04, 0.03, 0.02, 0. 
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However, it is more common in applications to fix parameters such as the damping coefficient v and vary the 
supercriticality parameter r. Then the results of Section V C 2 may be interpreted as indicating that, at least for small 
values of v, there are stable rolls provided 2r < u, i.e., close enough to onset. Thus the damped Nikolaevskiy equation 
differs fundamentally from the undamped version in that there is not a direct transition to soft-mode turbulence at 
onset when v > 0; instead, there is a sequence of bifurcations beginning with some initially stable roll state. Of 
course, when v is small this bifurcation sequence may occur very close to onset (i.e., for r — 0(v)). 

Our work complements a recent study by Tribelsky |12j . but differs in some significant respects. The first and 
most important difference is our discovery and analysis of the cusp of small-amplitude stable rolls, as described in 
Section V C 2 These rolls are the last to be destabilized as the damping is reduced, and so are particularly significant. 
Another difference is that our numerical study of the secondary stability problem has shown that rolls in the central 
wedge of apparent stability, which is predicted both by ourselves and by Tribelsky using asymptotic arguments, are in 
fact unstable to short-wave modes, except for exceedingly small values of v. So this central wedge (see Section VC 1 1 
in fact forms part of the stability boundary only for very small damping. 

Finally, we reiterate that our choice of damping term in ([6]) is by no means unique: indeed, Tribelsky |12j made 
a different choice, as discussed in Section |III| We conclude by noting that we have also computed numerically the 
secondary stability boundaries corresponding to those shown in Figure [2j but for Tribelsky's version of the damped 
Nikolaevskiy equation, and find similar results. 



[1] I. A. Beresnev and V. N. Nikolaevskiy, Physica D 66, 1-6 (1993). 

[2] V. N. Nikolaevskiy, Springer Lecture Notes in Engineering, 39, 210-221 (1989) 

[3] H. Fujisaka and T. Yamada, Prog. Theor. Phys. 106, 315-322 (2001). 

[4] D. Tanaka, Phys. Rev. E 70, 015202(R) (2004). 

[5] Y. Hidaka, J.-H. Huh, K.-i. Hayashi, S. Kai and M. I. Tribelsky, Phys. Rev. E 56, R6256-R6259 (1997). 

[6] M. I. Tribelsky and K. Tsuboi, Phys. Rev. Lett. 76, 1631-1634 (1996). 

[7] P. C. Matthews and S. M. Cox, Phys. Rev. E 62, R1473-R1476 (2000). 

[8] M. I. Tribelsky and M. G. Velarde, Phys. Rev. E 54, 4973-4981 (1992). 

[9] H. Sakaguchi and D. Tanaka, Phys. Rev. E 76, 025201(R) (2007). 
[10] H. Fujisaka, T. Honkawa and T. Yamada, Prog. Theor. Phys. 109, 911-918 (2003). 
[11] Y. Hidaka, K. Tamura, S. Kai, Prog. Theor. Phys. Supp. 161, 1-11 (2006). 
[12] M. I. Tribelsky, Preprint (2006). 

[13] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence. Springer, Berlin (1984). 
[14] G. I. Sivashinsky, 4 1177-1206 (1977). 

[15] T. Bansagi, D. Horvath, A. Toth, J. Yang, S. Kalliadasis and A. De Wit, Phys. Rev. E 68, 055301(R) (2003). 
[16] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851-1112 (1993). 
[17] D. Tanaka and Y. Kuramoto, Phys. Rev. E 68, 026219 (2003). 

[18] S. I. Anisimov, M. I. Tribel'skii, Ya. G. Epel'baum, Sov. Phys. JETP 51, 802-806 (1980). 

[19] J. D'Hernoncourt, S. Kalliadasis and A. De Wit, J. Chem. Phys. 123, 234503 (2005). 

[20] S. Kalliadasis, J. Yang and A. De Wit, Phys. Fluids 16, 1395-1409 (2004). 

[21] P. Brunet, Phys. Rev. E 76, 017204 (2007). 

[22] C. Misbah and A. Valance, Phys. Rev. E 49, 166-183 (1994). 

[23] A Doelman, PhD thesis, University of Utrecht (1990). 

[24] R. B. Hoyle, Phys. Rev. E 58, 7315-7318 (1998). 

[25] F. J. Mancebo and J. M. Vega, J. Fluid Mech. 560, 369-393 (2006). 

[26] S. M. Cox and P. C. Matthews, J. Comp. Phys. 176, 430-455 (2002). 



